rm(list=ls(all=TRUE))


library(doMPI)
library(MASS)
library(AER)
library(foreach)
library(spdep)
library(dplyr)
library(sphet)
library(ivpack)
set.seed(12345)
seeds = sample(seq(1:2000000),864000,replace  = F)


cl <- startMPIcluster()
registerDoMPI(cl)



spatialIV = function(N, beta, iv.coef, rho.x, rho.y, rho.z, cor,error, misspec){
    x.coord = runif(N,0,100)
    y.coord = runif(N,0,100)
    
    points = cbind(x.coord,y.coord)
#### create wrong W,
    pointsWrong <- cbind(runif(N,0,100),runif(N,0,100))


    dnb <- knearneigh(pointsWrong, 5, longlat = NULL)
    dsts <- knn2nb(dnb)
    lw <- nb2listw(dsts, style="B")
    wrongW <- listw2mat(lw)

    ## data generating W
    dnb <- knearneigh(points, 5, longlat = NULL)
    dsts <- knn2nb(dnb)
    lw <- nb2listw(dsts, style="B", zero.policy=TRUE)
    dgpW <- listw2mat(lw)

### create w for estimation
#### misspec is probability that cell value comes from DGP W
    multiplier <- matrix(rbinom(N *N, 1, misspec), ncol = N, nrow = N)

    estimationW <- matrix(NA, nrow = , N, ncol = N)
    estimationW[multiplier == 0] <- dgpW[multiplier == 0]
    estimationW[multiplier == 1] <- wrongW[multiplier == 1]

### row standardizing both matrices
    estimationW <- estimationW/rowSums(estimationW)
    dgpW <- dgpW/rowSums(dgpW)


    exogenous = matrix(rnorm(N*2), nrow = N, ncol = 2)

    inv_z = solve(diag(N) - (rho.z * (dgpW)))
    inv_x = solve(diag(N) - (rho.x * (dgpW)))
    inv_y = solve(diag(N) - (rho.y *(dgpW)))

    u = rnorm(N)

    instr = inv_z%*%u 
    gamma1 = c( 3, -2.5)
    gamma2 = c( -3 ,2.5)


    omega = matrix(c(1,cor,cor,1),nrow=2)
    tau = diag(c(error,error))
    covar = t(tau)%*%omega%*%tau  
    epsilon = mvrnorm(N, c(0,0), covar)

    x = inv_x%*%(( 2 + exogenous%*%gamma1 + instr*iv.coef + epsilon[,1]))
    y = inv_y%*%((-2 + exogenous%*%gamma2 + x*beta + epsilon[,2]))

                                        #  ## simple ols
    linear = lm(y ~ exogenous + x)  
    ## standard IV
    instrument = ivreg(y ~ exogenous + x | exogenous + instr)
    ## standard Iv with robust se
    iv.robustSE = robust.se(instrument)

    ## SAR-IV
    result_s2sls <- tryCatch({
        spiv = spreg(y ~  exogenous,listw = estimationW, endog =~ x, instruments = ~ instr, lag.instr = TRUE, initial.value = 0.2, model = "lag", het = TRUE)
        c( coefficients(spiv)[,1]["x"], sqrt(summary(spiv)$var[4,4]))
    },
    error=function(cond){
        result_s2sls <- c(NA, NA)
        return(result_s2sls)
    }
    )


### record correlation of spatial lag with correct and misspec W
    lagY <- dgpW %*%y

    lagYmis <- estimationW %*%y

### RECORD THIS
    correlation_spatLag <- cor(lagY, lagYmis)


### save coefficient estimates and SE for each method
    results = data.frame(c(summary(linear)$coefficients["x","Estimate"],summary(instrument)$coefficients["x","Estimate"], summary(instrument)$coefficients["x","Estimate"], result_s2sls[1]),  c(summary(linear)$coefficients["x","Std. Error"],summary(instrument)$coefficients["x","Std. Error"], iv.robustSE["x",2], result_s2sls[2]), c("ols", "2sls","2sls-robust",  "SAR-iv"), rep(summary(instrument,  diagnostics = TRUE)$diagnostics[1,3],4) , rep(rho.x,4), rep(rho.z,4), rep(rho.y, 4), rep(cor,4), rep(iv.coef,4), rep(beta,4), rep(error,4), rep(N,4), rep(misspec, 4), rep(correlation_spatLag, 4))
    return(results)
}




N = c(50,200)
cor = c(-0.5,0,0.5)
rhoy = seq(0,0.6, 0.3)
IV_strength = c(0.75, 1.5)
rhoz = seq(0,0.6, 0.3)
rhox = 0
error = 2#c(,3)
missp = c(0, 0.5, 1)


simulation_para = data.frame(do.call(rbind,replicate(1000,as.matrix(expand.grid(N,cor, IV_strength, rhoy,rhoz,rhox,error, missp)),simplify = F)))
names(simulation_para) = c("N","cor","IV.strength","rhoy", "rhoz", "rhox", "error", "missp")
simulation_para = simulation_para[order(simulation_para$N,simulation_para$cor, simulation_para$IV.strength, simulation_para$rhoy,simulation_para$rhoz,simulation_para$rhox, simulation_para$error, simulation_para$missp),]
dim(simulation_para)







writeLines(c(""), "log_test1.txt")

simulation_output_RR1 = foreach(i = 1:81000,
                                .combine = bind_rows,
                                .export = c("simulation_para","spatialIV","seeds"),
                                .verbose =FALSE,
                                .errorhandling = "remove",
                                .packages = c("AER", "ivpack","MASS","spdep","sphet")
                                ) %dopar% {
    library(AER)
    library(spdep)
    library(sphet)
    library(ivpack)
    library(MASS)
    set.seed(seeds[i])
    res = spatialIV(N=simulation_para[i,"N"],beta=2, iv.coef = simulation_para[i,"IV.strength"],  rho.x = simulation_para[i,"rhox"], rho.y = simulation_para[i,"rhoy"], rho.z = simulation_para[i,"rhoz"], cor = simulation_para[i,"cor"], error = simulation_para[i,"error"], misspec = simulation_para[i, "missp"])
    sink("log_test1.txt", append=TRUE)
    cat(paste("iter",i,"\n"))
    sink()
    return(res)
    rm(res)
                                    }
names(simulation_output_RR1) = c("Coef","SE","Method","f-stat", "rho.x", "rho.z", "rho.y", "cor","iv.coef", "treat", "Error", "N", "misspec", "correlation_spatial_lag")
save(simulation_output_RR1, file = "simulation_RR1.rda")

closeCluster(cl)
mpi.finalize()

